/* Mario Orsi - Soton Uni - 2005 
 definitions.h: constants, useful macros for plain, vector and math, and some file stuff.
 note: a matrix is treated as a linear array, e.g.:
 M = m[0] m[3] m[6]
     m[1] m[4] m[7]
     m[2] m[5] m[8]
*/

#ifndef DEFINITIONS_H
#define DEFINITIONS_H

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <signal.h>
#include <sys/time.h>


#define Sqr( x )     ((x) * (x))
#define Cube( x )    ((x) * (x) * (x))

/* Fundamental constants '_' denotes MKS units.  SHOULD NEVER BE ALTERED      */
#define	PI	3.14159265358979323846
#define DTOR	( PI / 180.0 )
#define FREE   3.0

/* Values from CODATA  */
#define	kB		1.3806505e-23 /* [J / K] */
#define kB_cal_K        3.29766597002e-24 /* [ cal / K ] */
#define AVOGAD		6.0221367e23 /* [mol^(-1)] */
#define kB_cal_K_mol    ( kB_cal_K * AVOGAD ) /* [ cal / (K * mol) ] */
#define AMU		1.6605402e-27 /* [kg] */
#define AMUSTR          "1.6605402e-27"
#define RTAMU           4.07497263794495e-14	/* [sqrt(amu)] */
#define RGAS		( AVOGAD * kB ) /* ...(8.314511)... [J / ( mol K )] */
#define EPS0            8.854187817e-12 /* [F/m] = [C^2/(Nm^2)] */
#define DEBYE_IN_Cm     ( 3.33564e-30 ) /* [Cm] - 1 Debye equals 3.33564*10^-30 coulomb*meter */
#define e_IN_COULOMB     ( 1.60217653e-19 ) /* [C] - "e" (charge of a proton) equals 1.602*10^-19 coulomb */


/* BRAHMS reduced units convertion factors - if multiplied by a reduced quantity they give that quantity in the unit they say */
// fundamental units:
#define MASS_amu                  18.01 /* to be read: MASS of water / amu = 18.01 */
#define MASS_kg                   2.987e-26 /* [kg] */
#define LENGTH_A                  3.035 /* [A] */
#define ENERGY_kcal_mol           0.152 /* [kcal/mol] */
// derived units:
#define LENGTH_m                  (LENGTH_A * 1e-10) /* [m] */
#define LENGTH_nm                 (LENGTH_A * 1e-1) /* [nm] */
#define VOLUME_nm3                (Cube( LENGTH_nm ))
#define VOLUME_m3                 (Cube( LENGTH_m ))
#define AREA_nm2                  (Sqr( LENGTH_nm ))
#define AREA_m2                   (Sqr( LENGTH_m ))
#define LENGTH_cm                 (LENGTH_m * 1e2)  /* [cm] */
#define ENERGY_cal_mol            ( ENERGY_kcal_mol * 1000 ) /* [cal/mol] */
// #define ENERGY_J                  1.056050421e-21  /* [J] */
#define J_per_cal                 ( 4.184 ) // J
#define ENERGY_J                  ( ENERGY_cal_mol * J_per_cal / AVOGAD ) 
#define FORCE_kcal_mol_nm        ( ENERGY_kcal_mol / LENGTH_nm )
#define FORCE_kJ_mol_nm          ( ENERGY_kcal_mol * J_per_cal / LENGTH_nm )
#define FORCE_N                   ( ENERGY_J /  LENGTH_m ) // Force [N]
#define FORCE_pN                  ( FORCE_N * 1e12 ) // Force [pN]
#define PRESSURE_atm              ( ENERGY_J / Cube( LENGTH_m ) / 101325 ) /* given that 1 atm = 101325 N/m^2 */
#define TEMPERATURE_K             ( ENERGY_J / kB )
#define TIME_s                    ( 1. / ( sqrt( ENERGY_J / ( MASS_kg * Sqr( LENGTH_m ) ) ) ) )
#define TIME_mus                  ( TIME_s * 1e6 ) // microsecond
#define TIME_ps                   ( TIME_s * 1e12 ) 
#define TIME_ns                   ( TIME_s * 1e9 ) 
#define TIME_fs                   ( TIME_s * 1e15 ) 
#define DENSITY_rel               ( MASS_kg / Cube( LENGTH_m ) / 1000 ) /* relative density = density in g/cm^3*/
#define DIPOLE_Cm                 ( sqrt( 4. * PI * EPS0 * Cube( LENGTH_m ) * ENERGY_J ) ) 
#define DIPOLE_D                  ( DIPOLE_Cm / DEBYE_IN_Cm ) 
#define CHARGE_C                  ( sqrt( 4. * PI * EPS0 * LENGTH_m * ENERGY_J ) )
#define RIGIDITY_kcal_molASq      ( ENERGY_kcal_mol / Sqr( LENGTH_A ) )  
#define TENSION_dyn_cm            ( FORCE_N * 1e3 / LENGTH_m ) // surface tension unit: dyn/cm = mN/m      - 1 dyn = 1e-5 N
//#define TENSION_dyn_cm            ( FORCE_N * 1e5 / LENGTH_CM ) // surface tension unit: dyn/cm = mN/m      - 1 dyn = 1e-5 N

// constants in reduced units
#define BETA               ( 4.6e-5 * PRESSURE_atm ) // isothermal compressibility of water, used in Berendsen's NPT technique 
                                                     // for lipid it should really be 4.5e-5, so its basically the same

/* miscellany */
#define T_room_K             293  // room temperature [K], same as 20 Celsius
#define kBT_room_J           ( kB * T_room_K ) // thermal energy at room temperature ~ 4.0453E-21 J
#define WATNUMDENSITY      0.03346 /* [molecules/A^2] - experimental number density of H20 */
#define HYDROTHICKNESS_DMPC    ( 25.4 / LENGTH_A ) // hydrocarbon chain region thickness - from Kucerka et al., Biophys J 
#define MAX_SITES_PER_LIPID    (50)  // Used in latDiff measurements, maybe it will need to be modified (increased) for very big lipids... 
#define OXYGEN_CHARGE_SSD       ( - 0.846 * e_IN_COULOMB ) // [C] effective charge computed from the SSD dipole value
#define HYDROGEN_CHARGE_SSD     ( - OXYGEN_CHARGE_SSD / 2. ) // [C]  
#define CHOL_CHARGE       (  1. * e_IN_COULOMB ) // [C] 
#define PHOS_CHARGE       ( -1. * e_IN_COULOMB ) // [C] 
#define GLY_PLUS_CHARGE       (  .2 * e_IN_COULOMB ) // [C] effective charges yielding a dipole of 4.2D when spaced 3A apart 
#define GLY_MINUS_CHARGE       (  -GLY_PLUS_CHARGE ) // [C] effective charges yielding a dipole of 4.2D when spaced 3A apart 
#define WAT_DIPOLE            ( 2.42 * DEBYE_IN_Cm ) // [C*m] from the ssd dipole which is 2.42 D
#define DMPC_BILAYER_BENDING_MODULUS_J   ( 0.56e-19 ) // [J] from Rawicz et al., Biophys J [2000]
#define DOPC_BILAYER_BENDING_MODULUS_J   ( 0.76e-19 ) // [J] from Pan et al., BJ [2008]
#define c0_DOPC_per_nm       ( -0.115 ) // [1/nm] DOPC spontaneous curvature, from Chen & Rand, BJ, [1997]
#define P1_DOPC_kBT_per_A    ( 0.5 * DOPC_BILAYER_BENDING_MODULUS_J / kBT_room_J * c0_DOPC_per_nm / 10 ) // pressure profile's first moment

/* text characters used somewhere around */
#define CHAR_MINUS  '-'
#define CHAR_ZERO   '0' /* denotes the character '0' in function SetupFiles */

#define	L_name		128			/* Max Length of file names  */
#define MAX( a, b )        ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) )  
#define MAXOF3( x, y, z )	MAX( x, MAX( y, z ) )
#define SUMSQ( x )	( x[ 0 ] * x[ 0 ] + x[ 1 ] * x[ 1 ] + x[ 2 ] * x[ 2 ] )
#define SUMSQ2( x )	( x[ 1 ] * x[ 1 ] + x[ 2 ] * x[ 2 ] + x[ 3 ] * x[ 3 ] )
#define	POW2( x )		( ( x ) * ( x ) )
#define POW3( x )		( ( x ) * ( x ) * ( x ) )
#define DIST( x, y )	sqrt( POW2( x[ 0 ] - y[ 0 ] ) + POW2( x[ 1 ] - y[ 1 ] ) + POW2( x[ 2 ] - y[ 2 ] ) )
#define	SIGN( x, y )	( ( y ) > 0 ? fabs( x ) : -fabs( x ) )
#define	ROUND( x )  ( int )( x + 0.5 )

#define Sgn( x, y )  ( ( ( y ) >= 0 ) ? ( x ) : ( - ( x ) ) )
#define IsEven( x )  ( ( x ) & ~1 )
#define IsOdd( x )   ( ( x ) & 1 )
#define Nint( x )                                             \
   ( ( ( x ) < 0. ) ? ( - ( int ) ( 0.5 - ( x ) ) ) : ( ( int ) ( 0.5 + ( x ) ) ) )
#define Min( x1, x2 )  ( ( ( x1 ) < ( x2 ) ) ? ( x1 ) : ( x2 ) )
#define Max(x1, x2)  (((x1) > (x2)) ? (x1) : (x2))
#define Min3(x1, x2, x3) \
   (((x1) < (x2)) ? (((x1) < (x3)) ? (x1) : (x3)) :         \
                    (((x2) < (x3)) ? (x2) : (x3)))
#define Max3(x1, x2, x3) \
   (((x1) > (x2)) ? (((x1) > (x3)) ? (x1) : (x3)) :         \
                    (((x2) > (x3)) ? (x2) : (x3)))
#define Clamp(x, lo, hi)                                    \
   (((x) >= (lo) && (x) <= (hi)) ? (x) :                    \
   (((x) < (lo)) ? (lo) : (hi)))

//typedef float real;
typedef double real;
typedef struct { real x, y, z; } VecR; /* vector with real components */
typedef struct { int x, y, z; } VecI;  /* vector with int components */
typedef struct { real x, y; } VecRTwo; /* two-dimensional vector with real components */
typedef struct { int i1, i2; } PairInt;  /* two-dimensional structure with two int components */


/* basic vector operations */
#define VSet(v, sx, sy, sz)                                 \
   (v).x = sx,                                              \
   (v).y = sy,                                              \
   (v).z = sz
#define VSetTwo(v, sx, sy)  (v).x = sx,  (v).y = sy
#define VCopy(v1, v2)                                       \
   (v1).x = (v2).x,                                         \
   (v1).y = (v2).y,                                         \
   (v1).z = (v2).z

#define VCopyTwo(v1, v2) (v1).x = (v2).x, (v1).y = (v2).y

#define VScale(v, s)                                        \
   (v).x *= s,                                              \
   (v).y *= s,                                              \
   (v).z *= s

#define VScaleTwo(v, s)   (v).x *= s, (v).y *= s

#define VSCopy(v2, s1, v1)                                  \
   (v2).x = (s1) * (v1).x,                                  \
   (v2).y = (s1) * (v1).y,                                  \
   (v2).z = (s1) * (v1).z

#define VSCopyTwo(v2, s1, v1)  (v2).x = (s1) * (v1).x, (v2).y = (s1) * (v1).y

#define VAdd(v1, v2, v3)                                    \
   (v1).x = (v2).x + (v3).x,                                \
   (v1).y = (v2).y + (v3).y,                                \
   (v1).z = (v2).z + (v3).z

#define VAddTwo(v1, v2, v3)   (v1).x = (v2).x + (v3).x, (v1).y = (v2).y + (v3).y

#define VSub(v1, v2, v3)                                    \
   (v1).x = (v2).x - (v3).x,                                \
   (v1).y = (v2).y - (v3).y,                                \
   (v1).z = (v2).z - (v3).z

#define VSubTwo(v1, v2, v3)   (v1).x = (v2).x - (v3).x, (v1).y = (v2).y - (v3).y

#define VMul(v1, v2, v3)                                    \
   (v1).x = (v2).x * (v3).x,                                \
   (v1).y = (v2).y * (v3).y,                                \
   (v1).z = (v2).z * (v3).z
#define VMulTwo(v1, v2, v3)                                    \
   (v1).x = (v2).x * (v3).x,                                \
   (v1).y = (v2).y * (v3).y
#define VDiv(v1, v2, v3)                                    \
   (v1).x = (v2).x / (v3).x,                                \
   (v1).y = (v2).y / (v3).y,                                \
   (v1).z = (v2).z / (v3).z
#define VDivTwo(v1, v2, v3)                                 \
   (v1).x = (v2).x / (v3).x,                                \
   (v1).y = (v2).y / (v3).y  

// Sum (square v1 and divide by v2, element by element) - useful to compute rotational kinEnergy 
#define VSumSqDiv(v1, v2)                                   \
   ((v1).x * (v1).x / (v2).x + (v1).y * (v1).y / (v2).y + (v1).z * (v1).z / (v2).z) 

#define VSAdd(v1, v2, s3, v3)                               \
   (v1).x = (v2).x + (s3) * (v3).x,                         \
   (v1).y = (v2).y + (s3) * (v3).y,                         \
   (v1).z = (v2).z + (s3) * (v3).z
#define VSAddTwo(v1, v2, s3, v3)                               \
   (v1).x = (v2).x + (s3) * (v3).x,                         \
   (v1).y = (v2).y + (s3) * (v3).y
#define VSSAdd(v1, s2, v2, s3, v3)                          \
   (v1).x = (s2) * (v2).x + (s3) * (v3).x,                  \
   (v1).y = (s2) * (v2).y + (s3) * (v3).y,                  \
   (v1).z = (s2) * (v2).z + (s3) * (v3).z
#define VSSSAdd(v1, s2, v2, s3, v3, s4, v4)                          \
   (v1).x = (s2) * (v2).x + (s3) * (v3).x + (s4) * (v4).x,                  \
   (v1).y = (s2) * (v2).y + (s3) * (v3).y + (s4) * (v4).y,                  \
   (v1).z = (s2) * (v2).z + (s3) * (v3).z + (s4) * (v4).z
#define VSSSSAdd(v1, s2, v2, s3, v3, s4, v4, s5, v5)                          \
   (v1).x = (s2) * (v2).x + (s3) * (v3).x + (s4) * (v4).x + (s5) * (v5).x,                  \
   (v1).y = (s2) * (v2).y + (s3) * (v3).y + (s4) * (v4).y + (s5) * (v5).y,                  \
   (v1).z = (s2) * (v2).z + (s3) * (v3).z + (s4) * (v4).z + (s5) * (v5).z

/* product operations */
#define VDot(v1, v2)                                        \
   ((v1).x * (v2).x + (v1).y * (v2).y + (v1).z * (v2).z)
#define VDotTwo( v1, v2 )  ( ( v1 ).x * ( v2 ).x + ( v1 ).y * ( v2 ).y )
#define VWDot(v1, v2, v3)                                   \
   ((v1).x * (v2).x * (v3).x + (v1).y * (v2).y * (v3).y +   \
   (v1).z * (v2).z * (v3).z)
#define VCross(v1, v2, v3)                                  \
   (v1).x = (v2).y * (v3).z - (v2).z * (v3).y,              \
   (v1).y = (v2).z * (v3).x - (v2).x * (v3).z,              \
   (v1).z = (v2).x * (v3).y - (v2).y * (v3).x

/* matrix-vector products (for both the matrix and its transpose) */
#define MVMul(v1, m, v2)                                         \
   (v1).x = (m)[0] * (v2).x + (m)[3] * (v2).y + (m)[6] * (v2).z, \
   (v1).y = (m)[1] * (v2).x + (m)[4] * (v2).y + (m)[7] * (v2).z, \
   (v1).z = (m)[2] * (v2).x + (m)[5] * (v2).y + (m)[8] * (v2).z
#define MVMulT(v1, m, v2)                                        \
   (v1).x = (m)[0] * (v2).x + (m)[1] * (v2).y + (m)[2] * (v2).z, \
   (v1).y = (m)[3] * (v2).x + (m)[4] * (v2).y + (m)[5] * (v2).z, \
   (v1).z = (m)[6] * (v2).x + (m)[7] * (v2).y + (m)[8] * (v2).z

/* other useful operations */
#define MCopy( m1, m2 )   \
(m1)[0] = (m2)[0], (m1)[1] = (m2)[1], (m1)[2] = (m2)[2], \
(m1)[3] = (m2)[3], (m1)[4] = (m2)[4], (m1)[5] = (m2)[5], \
(m1)[6] = (m2)[6], (m1)[7] = (m2)[7], (m1)[8] = (m2)[8]

#define MScale(m, s)                                       \
   (m)[0] *= s,(m)[3] *= s,(m)[6] *= s,                    \
   (m)[1] *= s,(m)[4] *= s,(m)[7] *= s,                    \
   (m)[2] *= s,(m)[5] *= s,(m)[8] *= s

/* other useful operations */
#define MExpScale(m1, s, m2)                                       \
   (m1)[0] *= exp(s * (m2)[0]), (m1)[3] *= exp(s * (m2)[3]), (m1)[6] *= exp(s * (m2)[6]),                    \
   (m1)[1] *= exp(s * (m2)[1]), (m1)[4] *= exp(s * (m2)[4]), (m1)[7] *= exp(s * (m2)[7]),                    \
   (m1)[2] *= exp(s * (m2)[2]), (m1)[5] *= exp(s * (m2)[5]), (m1)[8] *= exp(s * (m2)[8])

// Matrix Determinant, assuming a 3 x 3 matrix
#define MatDet( m )                                         \
   (m)[0]*(m)[4]*(m)[8] - (m)[0]*(m)[7]*(m)[5] - (m)[3]*(m)[1]*(m)[8] + (m)[3]*(m)[7]*(m)[2] + (m)[6]*(m)[1]*(m)[5] - (m)[6]*(m)[4]*(m)[2]  

// Matrix addition - element by element
#define MAdd( m1, m2, m3 ) \
(m1)[0] = (m2)[0] + (m3)[0],  \
(m1)[1] = (m2)[1] + (m3)[1],  \
(m1)[2] = (m2)[2] + (m3)[2],  \
(m1)[3] = (m2)[3] + (m3)[3],  \
(m1)[4] = (m2)[4] + (m3)[4],  \
(m1)[5] = (m2)[5] + (m3)[5],  \
(m1)[6] = (m2)[6] + (m3)[6],  \
(m1)[7] = (m2)[7] + (m3)[7],  \
(m1)[8] = (m2)[8] + (m3)[8]

// Matrix subtraction - element by element
#define MSub( m1, m2, m3 ) \
(m1)[0] = (m2)[0] - (m3)[0],  \
(m1)[1] = (m2)[1] - (m3)[1],  \
(m1)[2] = (m2)[2] - (m3)[2],  \
(m1)[3] = (m2)[3] - (m3)[3],  \
(m1)[4] = (m2)[4] - (m3)[4],  \
(m1)[5] = (m2)[5] - (m3)[5],  \
(m1)[6] = (m2)[6] - (m3)[6],  \
(m1)[7] = (m2)[7] - (m3)[7],  \
(m1)[8] = (m2)[8] - (m3)[8]

// Matrix accumulation sum: m1 += m2 (element by element)
#define MMAdd( m1, m2 )   \
(m1)[0] += (m2)[0],   \
(m1)[1] += (m2)[1],   \
(m1)[2] += (m2)[2],   \
(m1)[3] += (m2)[3],   \
(m1)[4] += (m2)[4],   \
(m1)[5] += (m2)[5],   \
(m1)[6] += (m2)[6],   \
(m1)[7] += (m2)[7],   \
(m1)[8] += (m2)[8]

// Matrix accumulation subtraction: m1 -= m2 (element by element)
#define MMSub( m1, m2 )   \
(m1)[0] -= (m2)[0],   \
(m1)[1] -= (m2)[1],   \
(m1)[2] -= (m2)[2],   \
(m1)[3] -= (m2)[3],   \
(m1)[4] -= (m2)[4],   \
(m1)[5] -= (m2)[5],   \
(m1)[6] -= (m2)[6],   \
(m1)[7] -= (m2)[7],   \
(m1)[8] -= (m2)[8]

// Matrix addition with scaling - element by element
#define MSAdd( m1, m2, s3, m3 ) \
(m1)[0] = (m2)[0] + (s3)*(m3)[0],  (m1)[3] = (m2)[3] + (s3)*(m3)[3],  (m1)[6] = (m2)[6] + (s3)*(m3)[6],  \
(m1)[1] = (m2)[1] + (s3)*(m3)[1],  (m1)[4] = (m2)[4] + (s3)*(m3)[4],  (m1)[7] = (m2)[7] + (s3)*(m3)[7],  \
(m1)[2] = (m2)[2] + (s3)*(m3)[2],  (m1)[5] = (m2)[5] + (s3)*(m3)[5],  (m1)[8] = (m2)[8] + (s3)*(m3)[8]

#define MMSAdd( m1, s2, m2 )  MSAdd( m1, m1, s2, m2 )

#define MZero( m ) (m)[0] = 0, (m)[1] = 0, (m)[2] = 0, (m)[3] = 0, (m)[4] = 0, (m)[5] = 0, (m)[6] = 0, (m)[7] = 0, (m)[8] = 0

// Set off-diagonal elements to zero
#define MZeroOffDiags( m )   (m)[1] = 0, (m)[2] = 0, (m)[3] = 0, (m)[5] = 0, (m)[6] = 0, (m)[7] = 0

// Trace operator
#define Trace( m ) ( (m)[0] +  (m)[4] +  (m)[8] )

// macro to extract (pL - pN) from the stress/pressure tensor
#define ExtractLatPress( m ) ( 0.5 * ( (m)[0] +  (m)[4] ) -  (m)[8] )

// macro to extract (pN - pL) from the stress/pressure tensor
#define ExtractNormPress( m ) ( (m)[8] - 0.5 * ( (m)[0] +  (m)[4] )  )

// macro to extract (pL - pN) from the stress/pressure tensor
#define ExtractLatPressComponent( m ) ( 0.5 * ( (m)[0] +  (m)[4] ) )

// macro to extract (pN - pL) from the stress/pressure tensor
#define ExtractNormPressComponent( m ) ( (m)[8] )


// Outer product between vectors v1, v2 - result is (stored in) matrix m // checked 12-09-2006
#define VOuter( m, v1, v2 )   \
(m)[0] = (v1).x * (v2).x,     \
(m)[1] = (v1).y * (v2).x,     \
(m)[2] = (v1).z * (v2).x,     \
(m)[3] = (v1).x * (v2).y,     \
(m)[4] = (v1).y * (v2).y,     \
(m)[5] = (v1).z * (v2).y,     \
(m)[6] = (v1).x * (v2).z,     \
(m)[7] = (v1).y * (v2).z,     \
(m)[8] = (v1).z * (v2).z  

#define VProd(v)                                            \
   ((v).x * (v).y * (v).z)
#define VGe(v1, v2)                                         \
   ((v1).x >= (v2).x && (v1).y >= (v2).y && (v1).z >= (v2).z)
#define VLt(v1, v2)                                         \
   ((v1).x < (v2).x && (v1).y < (v2).y && (v1).z < (v2).z)
#define VLinear(p, s)                                       \
   (((p).z * (s).y + (p).y) * (s).x + (p).x)
#define VSetAll(v, s)                                       \
   VSet (v, s, s, s)
#define VSetAllTwo(v, s)   VSetTwo (v, s, s)
#define VAddCon(v1, v2, s)                                  \
   (v1).x = (v2).x + (s),                                   \
   (v1).y = (v2).y + (s),                                   \
   (v1).z = (v2).z + (s)
#define VCSum(v)                                            \
   ((v).x + (v).y + (v).z)
#define VComp(v, k)                                         \
   *((k == 0) ? &(v).x : ((k == 1) ? &(v).y : &(v).z))

/* two operations for converting between vectors and array elements */
#define VToLin(a, n, v)                                     \
   a[(n) + 0] = (v).x,                                      \
   a[(n) + 1] = (v).y,                                      \
   a[(n) + 2] = (v).z
#define VFromLin(v, a, n)                                   \
   VSet (v, a[(n) + 0], a[(n) + 1], a[(n) + 2])


/* additional vector operations */
#define VZeroTwo(v)  VSetAllTwo (v, 0)
#define VZero(v)  VSetAll (v, 0)
#define VLenSq(v)  VDot (v, v)
#define VLenSqTwo(v)  VDotTwo( v, v )
#define VWLenSq(v1, v2)  VWDot(v1, v2, v2)
#define VLen(v)  sqrt (VDot (v, v))
#define VVAdd(v1, v2)  VAdd (v1, v1, v2)
#define VVAddTwo(v1, v2)  VAddTwo (v1, v1, v2)
#define VVSub(v1, v2)  VSub (v1, v1, v2)
#define VVSAdd(v1, s2, v2) VSAdd (v1, v1, s2, v2)
#define VVSAddTwo(v1, s2, v2) VSAddTwo (v1, v1, s2, v2)
#define VInterp(v1, s2, v2, v3)                             \
   VSSAdd (v1, s2, v2, 1. - (s2), v3)

/* quaternions: definition of dataStructs and operations */
typedef struct {
   real u1, u2, u3, u4;
} Quat;

#define QSet(q, s1, s2, s3, s4)                             \
   (q).u1 = s1,                                             \
   (q).u2 = s2,                                             \
   (q).u3 = s3,                                             \
   (q).u4 = s4
#define QZero(q)  QSet (q, 0, 0, 0, 0)
#define QScale(q, s)                                        \
   (q).u1 *= s,                                             \
   (q).u2 *= s,                                             \
   (q).u3 *= s,                                             \
   (q).u4 *= s
#define QSAdd(q1, q2, s3, q3)                               \
   (q1).u1 = (q2).u1 + (s3) * (q3).u1,                      \
   (q1).u2 = (q2).u2 + (s3) * (q3).u2,                      \
   (q1).u3 = (q2).u3 + (s3) * (q3).u3,                      \
   (q1).u4 = (q2).u4 + (s3) * (q3).u4
#define QLenSq(q)                                           \
   (Sqr ((q).u1) + Sqr ((q).u2) + Sqr ((q).u3) +            \
   Sqr ((q).u4))
#define QMul(q1, q2, q3)                                    \
   (q1).u1 =   (q2).u4 * (q3).u1 - (q2).u3 * (q3).u2 +      \
               (q2).u2 * (q3).u3 + (q2).u1 * (q3).u4,       \
   (q1).u2 =   (q2).u3 * (q3).u1 + (q2).u4 * (q3).u2 -      \
               (q2).u1 * (q3).u3 + (q2).u2 * (q3).u4,       \
   (q1).u3 = - (q2).u2 * (q3).u1 + (q2).u1 * (q3).u2 +      \
               (q2).u4 * (q3).u3 + (q2).u3 * (q3).u4,       \
   (q1).u4 = - (q2).u1 * (q3).u1 - (q2).u2 * (q3).u2 -      \
               (q2).u3 * (q3).u3 + (q2).u4 * (q3).u4

/* complex variables */
typedef struct {
  real R, I;
} Cmplx;

#define CSet(a, x, y)                                       \
   a.R = x,                                                 \
   a.I = y
#define CAdd(a, b, c)                                       \
   a.R = b.R + c.R,                                         \
   a.I = b.I + c.I
#define CSub(a, b, c)                                       \
   a.R = b.R - c.R,                                         \
   a.I = b.I - c.I
#define CMul(a, b, c)                                       \
  a.R = b.R * c.R - b.I * c.I,                              \
  a.I = b.R * c.I + b.I * c.R

typedef struct {
  real u[ 9 ];
} RMat;

/* mar01
 constant used in PutConfig() note: 2^(15)=32768... */

#define SCALE_FAC 32767

/* feb22 */
#define PI 3.14159265358979323846

/* macro introduced to aid readability by hiding the fact that matrices
are in fact stored as singly-indexed array */
#define MAT( a, n, i, j )  ( a )[ ( i ) + n * ( j ) ]

/* macros handling allocations of one- and two-dimensional arrays of
   any kind of variable or structure */
#define AllocMem(a, n, t)  a = (t *) malloc ((n) * sizeof (t))
#define AllocMem2(a, n1, n2, t)                             \
   AllocMem (a, n1, t *);                                   \
   AllocMem (a[0], (n1) * (n2), t);                         \
   for (k = 1; k < n1; k ++) a[k] = a[k - 1] + n2;

#define MAX_MPEX_ORD  2
#define I(i, j)  ((i) * ((i) + 1) / 2 + (j))
#define c(i, j)  c[I(i, j)]
#define s(i, j)  s[I(i, j)]

typedef struct {
  real c[I(MAX_MPEX_ORD, MAX_MPEX_ORD) + 1], s[I(MAX_MPEX_ORD, MAX_MPEX_ORD) + 1];
} MpTerms;
typedef struct {
  MpTerms le, me;
  int occ;
} MpCell;

/* shorthand for frequently appearing loops */
#define DO_SITE          for ( n = 0; n < nSites; n++ )
#define DO_ATOM          for ( a = 0; a < nAtoms; a++ )
#define DO_WATER         for ( n = 0; n < nWaters; n++ )
#define DO_LIPID         for ( n = 0; n < nLipids; n++ )
#define DO_DMPCLIPID      for ( s = 0; s < dmpcStruct.nBeads; s++ ) 
#define DO_CELL( j, m )  for ( j = cellList[ m ]; j >= 0; j = cellList[ j ] )

/* macros dealing with Periodic Boundary Conditions - TO BE MODIFIED FOR NON-ISOSHAPE NPT DYNAMICS!!! */
#define VWrap(v, t)                                         \
   if (v.t >= 0.5 * region.t)      v.t -= region.t;         \
   else if (v.t < -0.5 * region.t) v.t += region.t

#define VWrapS(v, t)                                         \
   if (v.t >= 0.5 )      v.t -= 1.;         \
   else if (v.t < -0.5 ) v.t += 1.

#define VShift(v, t)                                        \
   if (v.t >= 0.5 * region.t)      shift.t -= region.t;     \
   else if (v.t < -0.5 * region.t) shift.t += region.t

#define VShiftWrap(v, t)                                    \
   if (v.t >= 0.5 * region.t) {                             \
     shift.t -= region.t;                                   \
     v.t -= region.t;                                       \
   } else if (v.t < -0.5 * region.t) {                      \
     shift.t += region.t;                                   \
     v.t += region.t;                                       \
   }

#define VCellWrap(t)                                        \
   if (m2v.t >= cells.t) {                                  \
     m2v.t = 0;                                             \
     shift.t = region.t;                                    \
   } else if (m2v.t < 0) {                                  \
     m2v.t = cells.t - 1;                                   \
     shift.t = - region.t;                                  \
   }

#define VWrapAll(v)                                         \
   {VWrap (v, x);                                           \
   VWrap (v, y);                                            \
   VWrap (v, z);}

#define VWrapXY( v ) { VWrap( v, x );  VWrap( v, y ); }

#define VWrapAllS(v)                                         \
   {VWrapS (v, x);                                           \
   VWrapS (v, y);                                            \
   VWrapS (v, z);}
#define VShiftAll(v)                                        \
   {VShift (v, x);                                          \
   VShift (v, y);                                           \
   VShift (v, z);}
#define VCellWrapAll()                                      \
   {VCellWrap (x);                                          \
   VCellWrap (y);                                           \
   VCellWrap (z);}

/* cell offsets used for interaction and neighbour-list computations */
#define OFFSET_VALS                                           \
   { {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0}, {-1,1,0},            \
     {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}, {-1,1,1}, {-1,0,1},  \
     {-1,-1,1}, {0,-1,1}, {1,-1,1}                            \
   }

#define N_OFFSET  14

/* measurements stuff */
typedef struct { /* property measurements */
  real val, /* actual measured value */ 
    sum,    /* sum accumulated, when averaged gives the mean */
    sum2;   /* sum of squares, when averaged gives the standard deviation */
} Prop;

#define PropZero(v)  v.sum = v.sum2 = 0. // initialisation
#define PropAccum(v)  v.sum += v.val, v.sum2 += Sqr (v.val)

// v.sum := mean, v.sum2 := sqrt( variance ) = standard deviation = rms fluctuation
#define PropAvg(v, n)  v.sum /= n, v.sum2 = sqrt (Max (v.sum2 / n - Sqr (v.sum), 0.))
#define PropEst(v)  v.sum, v.sum2

typedef struct {
  real time;
  int left, right, up, circAL, circAR, circBL, circBR, idA, idB;
} EvTree;

#define SITE_LIMIT  10000000

#define NameVal(x)                                          \
   if (! strncmp (bp, #x, strlen (#x))) {                   \
     bp += strlen (#x);                                     \
     x = strtod (bp, &bp);                                  \
   }

#define NHIST  (5)  /* NDIM + 2 */

/* file stuff */
#define ReadF( x )       fread( &x, sizeof( x ), 1, fp )
#define WriteF( x )      fwrite( &x, sizeof( x ), 1, fp )
#define ReadFN( x, n )   fread( x, sizeof( x[ 0 ] ), n, fp )
#define WriteFN( x, n )  fwrite( x, sizeof( x[ 0 ] ), n, fp )

#endif
